
# set working directory
PATH <- "~/Desktop/Research&Politics/replication-files"
# required libraries
library("tidyverse")
library("haven")
library("lme4")
library("gridExtra")
# utility functions
std <- function (x,select=TRUE) (x - mean(x[select],na.rm=TRUE)) / sd(x[select],na.rm=TRUE)
sim_data_glm <- function(m, x, effect_size) {
    invlogit <- function(x) exp(x)/(1+exp(x))
    res <- m@frame
    J <- fitted(m) %>% length()
    K <- ranef(m)[[1]] %>% nrow()

    sigma.a <- VarCorr(m)[[1]][1] %>% sqrt()
    mu      <- fixef(m)["(Intercept)"]
    b       <- fixef(m)
    # set effect size
    b[x] <- effect_size
    # create data
    raneff <- rnorm(K, 0, sigma.a)[as.integer(res$canton)]

    res[[names(res)[1]]] <- ((model.matrix(m) %*% b)[,1] + raneff) %>%
        invlogit() %>%
        rbinom(length(.), 1, .)
    return(res)
}
power_glm <- function(m, x = NULL, effect_size = NULL, n.sims = 1000) {
    signif <- rep(NA, n.sims)
    for (s in 1:n.sims) {
        print(s)
        fake      <- sim_data_glm(m, x = x, effect_size = effect_size)
        lme_power <- glmer(formula(m), data = fake, family = binomial(logit), nAGQ = 0)
        p_val     <- summary(lme_power)$coefficients[2, 4] 
        signif[s] <- p_val < 0.05
    }
    power <- mean(signif)
    return(power)
    # The function returns the proportion of the 1000 simulations where the result is statistically signiﬁcant; thus, the power (as computed via simulation) for a study with J children measured at K equally spaced times during the year.
}

## open data
d <- read_stata(file.path(PATH, "swiss-mps.dta"))
d <- d %>%
    mutate(
        index_direct_democracy_1_std = std(index_direct_democracy_1),
        index_direct_democracy_2_std = std(index_direct_democracy_2),
        pop_std                      = std(pop),
        unemployment_std             = std(unemployment),
        sender_occupation            = factor(sender_occupation),
        sender_name                  = factor(sender_name),
        batch                        = factor(batch),
        party                        = factor(party)
    )

# model specification
f <- list(
    replied ~ index_direct_democracy_1_std + language_french +
        female + sender_occupation + sender_union + sender_name + 
        batch + party +
        pop_std + unemployment_std + (1 | canton),
    replied ~ index_direct_democracy_2_std + language_french +
        female + sender_occupation + sender_union + sender_name + 
        batch + party +
        pop_std + unemployment_std + (1 | canton)
)

# power calculation for index 1 and 2
# WARNING: This is going to take some time...
for (i in 1:2) {
    m <- glmer(f[i], data = d, family = binomial(logit))

    set.seed(123456)
    fixef(m)[2]
    x <- seq(-0.3, 0.3, 0.05)
    x <- seq(-0.5, 0.5, 0.01)
    x <- seq(-0.5, 0.5, 0.02)
    p <- x %>%
        map_dbl(~power_glm(m, x = sprintf("index_direct_democracy_%s_std", i), effect_size = .x, n.sims = 1000))
    # save results from simulation
    write_rds(p, file.path(PATH, sprintf('res-%s.rds', i)))    
}

# create figure
p1 <- read_rds(file.path(PATH, "res-1.rds"))
p2 <- read_rds(file.path(PATH, "res-2.rds"))
x <- seq(-0.3, 0.3, 0.01)
x <- seq(-0.5, 0.5, 0.01)
x <- seq(-0.5, 0.5, 0.02)

g1 <- qplot(x, p1, geom = c("line", "point")) +
    scale_x_continuous("Effect Size") +
    scale_y_continuous("Statistical Power", limits = c(0, 1)) +
    labs(title = "Direct Democracy: Formal rules") +
    theme_bw() + theme(text = element_text(size = 8.5))
g2 <- qplot(x, p2, geom = c("line", "point")) +
    scale_x_continuous("Effect Size") +
    scale_y_continuous("Statistical Power", limits = c(0, 1)) +
    labs(title = "Direct Democracy: Actual use") +
    theme_bw() + theme(text = element_text(size = 8.5))
dev.new(width = 7.328859, height = 3.288591)
grid.arrange(g1, g2, ncol = 2)

# bivariate model
f <- replied ~ index_direct_democracy_1_std + (1 | canton)
m <- glmer(f, data = d, family = binomial(logit))
fixef(m)[2]
x <- seq(0, 0.3, 0.01)
p <- x %>%
    map_dbl(~power_glm(m, x = "index_direct_democracy_1_std", effect_size = .x, n.sims = 100))

qplot(x, p, geom = c("line", "point")) +
    scale_x_continuous("Effect Size") +
    scale_y_continuous("Power", limits = c(0, 1)) +
    theme_bw()
